#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<math.h>
#include<time.h>
#include<direct.h>
#include<sys\\stat.h>
#include<sys\\types.h>

#define mSol 1
#define mTerra 3.003167261575587e-6
#define mJupiter 9.543512141169373e-4

#define PI 3.141592653589793
#define G 2.959572553206070e-004
#define perielioTerra 0.983235919814466
#define vPerielioTerra 6.389661159796174
#define excentricidadeTerra 0.0167

void inicializacao(void);
void artigo2(void);
void oito(void);
void artigo1(void);
void coordenadasCentroMassa(void);
void criarArquivos(void);
void F3(double *dx, double *z, double *distancia, double *m, int dimensao);
void calcularDistancia(double *z,double *distancia, int dimensao, int corpos);
void jacobiana3C(double* j, double m1,double m2,double m3,double r1x,double r1y,double r2x,double r2y,double r3x,double r3y,int tamanho);
void copiarVetor(double* v,double* v2);
void produtoEscalarVetor(double* v,double a,int tamanho, double *resultado);
void produtoEscalar(double* v,double* v2,int tamanho, double* resultado);
void diferencaVetores(double* v,double* v2,int tamanho, double *resultado);
void somaVetores(double* v,double* v2,int tamanho, double *resultado);
void produtoMatrizVetor(double *J,double* delta,int tamanho, double *resultado);
void escreverArquivos(void);
void imprimir(double* v,int tamanho);
double calculaNorma(double *v,int dimensao);
void criarMatrizRK45(void);
void finalizar(void);

double *K1,*K2,*K3,*K4,*K5,*K6;
double  e,a2,*m,tFinal,dt,*a,*b4,*b5,*c2,*y,*Y,*norma,*temp,*temp2,*temp3,*temp4,t;
double  *dx, *distancia,tempo,*J,*z,*z1,*z2,tempo,tol=10e-8,delta;
int dimensao, corpos,escrever,tamanho,tamanho2,renormalizacao,orbitas=0,stages=6,flag=1;
char *exemplo,*aux,*SO;
clock_t inicio,fim,total;
FILE *f[8];

int main(int argc, char *argv[])
{
	double temp3;
	inicio=clock();

	SO=argv[1];
	escrever=atoi(argv[2]);
	e=atof(argv[3]);
	a2=atof(argv[4]);
	tFinal=atof(argv[5]);
	dt=atof(argv[6]);

	inicializacao();
	criarArquivos();
	escreverArquivos();
	criarMatrizRK45();

	//RK4

	while(flag)
	{
		//temp2=soma parcial

		//RK5
		//stage 1
		F3(temp, z, distancia, m, dimensao);
		produtoEscalarVetor(temp,dt,tamanho,K1);

		//stage 2
		produtoEscalarVetor(K1,a[1*stages+0],tamanho,temp);
		somaVetores(temp,z,tamanho,temp2);
		F3(temp, temp2, distancia, m, dimensao);
		produtoEscalarVetor(temp,dt,tamanho,K2);

		//stage 3
		produtoEscalarVetor(K1,a[2*stages+0],tamanho,temp);
		somaVetores(temp,z,tamanho,temp2);
		produtoEscalarVetor(K2,a[2*stages+1],tamanho,temp);
		somaVetores(temp,temp2,tamanho,temp2);
		F3(temp, temp2, distancia, m, dimensao);
		produtoEscalarVetor(temp,dt,tamanho,K3);

		//stage 4
		produtoEscalarVetor(K1,a[3*stages+0],tamanho,temp);
		somaVetores(temp,z,tamanho,temp2);
		produtoEscalarVetor(K2,a[3*stages+1],tamanho,temp);
		somaVetores(temp,temp2,tamanho,temp2);
		produtoEscalarVetor(K3,a[3*stages+2],tamanho,temp);
		somaVetores(temp,temp2,tamanho,temp2);
		F3(temp, temp2, distancia, m, dimensao);
		produtoEscalarVetor(temp,dt,tamanho,K4);

		//stage 5
		produtoEscalarVetor(K1,a[4*stages+0],tamanho,temp);
		somaVetores(temp,z,tamanho,temp2);
		produtoEscalarVetor(K2,a[4*stages+1],tamanho,temp);
		somaVetores(temp,temp2,tamanho,temp2);
		produtoEscalarVetor(K3,a[4*stages+2],tamanho,temp);
		somaVetores(temp,temp2,tamanho,temp2);
		produtoEscalarVetor(K4,a[4*stages+3],tamanho,temp);
		somaVetores(temp,temp2,tamanho,temp2);
		F3(temp, temp2, distancia, m, dimensao);
		produtoEscalarVetor(temp,dt,tamanho,K5);

		//stage 6
		produtoEscalarVetor(K1,a[5*stages+0],tamanho,temp);
		somaVetores(temp,z,tamanho,temp2);
		produtoEscalarVetor(K2,a[5*stages+1],tamanho,temp);
		somaVetores(temp,temp2,tamanho,temp2);
		produtoEscalarVetor(K3,a[5*stages+2],tamanho,temp);
		somaVetores(temp,temp2,tamanho,temp2);
		produtoEscalarVetor(K4,a[5*stages+3],tamanho,temp);
		somaVetores(temp,temp2,tamanho,temp2);
		produtoEscalarVetor(K5,a[5*stages+4],tamanho,temp);
		somaVetores(temp,temp2,tamanho,temp2);
		F3(temp, temp2, distancia, m, dimensao);
		produtoEscalarVetor(temp,dt,tamanho,K6);

		produtoEscalarVetor(K1,b4[0],tamanho,temp);
		somaVetores(z,temp,tamanho,z1);
		produtoEscalarVetor(K2,b4[1],tamanho,temp);
		somaVetores(z1,temp,tamanho,z1);
		produtoEscalarVetor(K3,b4[2],tamanho,temp);
		somaVetores(z1,temp,tamanho,z1);
		produtoEscalarVetor(K4,b4[3],tamanho,temp);
		somaVetores(z1,temp,tamanho,z1);
		produtoEscalarVetor(K5,b4[4],tamanho,temp);
		somaVetores(z1,temp,tamanho,z1);
		produtoEscalarVetor(K6,b4[5],tamanho,temp);
		somaVetores(z1,temp,tamanho,z1);

		produtoEscalarVetor(K1,b5[0],tamanho,temp);
		somaVetores(z,temp,tamanho,z2);
		produtoEscalarVetor(K2,b5[1],tamanho,temp);
		somaVetores(z2,temp,tamanho,z2);
		produtoEscalarVetor(K3,b5[2],tamanho,temp);
		somaVetores(z2,temp,tamanho,z2);
		produtoEscalarVetor(K4,b5[3],tamanho,temp);
		somaVetores(z2,temp,tamanho,z2);
		produtoEscalarVetor(K5,b5[4],tamanho,temp);
		somaVetores(z2,temp,tamanho,z2);
		produtoEscalarVetor(K6,b5[5],tamanho,temp);
		somaVetores(z2,temp,tamanho,z2);

		diferencaVetores(z1,z2,tamanho,temp);
		temp3=calculaNorma(temp,tamanho);
		temp3=temp3/dt;
		if(temp3<=tol)
		{
			t+=dt;
			copiarVetor(z,z1);
			escreverArquivos();
		}

		delta=tol/temp3;
		delta=pow(delta,0.25);
		delta=0.84*delta;

		/*if(delta<0.1)
		{
			dt=0.1*dt;
		}
		else*/ if(delta>4)
			 {
				dt=4*dt;
			 }
		     else
		     {
		    	 dt=delta*dt;
		     }

		if(t>=tFinal)
		{
			flag=0;
		}
		else if(t+dt>tFinal)
		{
			dt=tFinal-t;
		}
		//printf("%e\n",temp3);
	}

	int minutos=0;
	int segundos=0;
	int horas=0;
	tempo=0;
	int tempo2=0;

	fim=clock();
	total=fim-inicio;
	tempo=(double)total/CLOCKS_PER_SEC;
	tempo2=(int)tempo;
	horas=(tempo2/3600);
	minutos=(tempo2/60)-horas*60;
	segundos=(tempo2 % 60);
	printf("%d h %d min %d s\n",horas, minutos ,segundos);
	printf("%lf \n",tempo);

	sprintf(aux,"%0.10lf ",tempo);
	fputs(aux,f[2*corpos]);
	fputs("\n",f[2*corpos]);

	finalizar();

	return 0;
}

void inicializacao(void)
{
	t=0;
	dimensao=2;
	corpos=3;
	renormalizacao=1;
	tamanho=4*corpos*dimensao+2;
	tamanho2=2*corpos*dimensao;

	int i;

	m =(double *)calloc(corpos, sizeof(double));
	z =(double *)calloc(tamanho, sizeof(double));
	z1 =(double *)calloc(tamanho, sizeof(double));
	z2 =(double *)calloc(tamanho, sizeof(double));
	temp =(double *)calloc(tamanho, sizeof(double));
	temp2 =(double *)calloc(tamanho, sizeof(double));
	temp3 =(double *)calloc(1, sizeof(double));
	temp4 =(double *)calloc(1, sizeof(double));
	distancia =(double *)calloc(corpos*corpos, sizeof(double));
	J=(double *)calloc(pow(tamanho2,2), sizeof(double));
	K1=(double*)calloc(tamanho,sizeof(double));
	K2=(double*)calloc(tamanho,sizeof(double));
	K3=(double*)calloc(tamanho,sizeof(double));
	K4=(double*)calloc(tamanho,sizeof(double));
	K5=(double*)calloc(tamanho,sizeof(double));
	K6=(double*)calloc(tamanho,sizeof(double));
	aux =(char *)calloc(1000, sizeof(char));

	//delta inicial aleatorio
	srand( (unsigned)time(NULL) );
	for(i=tamanho2;i<tamanho-2;i++)
	{
		//delta[i]=rand()/(double)RAND_MAX;
		//delta[i]=delta[i]-0.5;
		z[i]=1.0;
	}

	// norma delta inicial nao normalizado
	*temp3=calculaNorma(z+tamanho2,tamanho2);
	//normalizacao delta inicial
	for(i=tamanho2;i<tamanho-2;i++)
	{
		z[i]=z[i]/temp3[0];
	}
	imprimir(z+tamanho2,tamanho2);

	artigo2();
	coordenadasCentroMassa();

	//velocidade centro de massa ser nula
	for(i=0;i<dimensao;i++)
	{
		z[corpos*dimensao+i]=-z[(corpos+1)*dimensao+i]-z[(corpos+2)*dimensao+i];
	}
}

void artigo2(void)
{
	exemplo="C:\\dados";
	mkdir(exemplo);
	exemplo="C:\\dados\\artigo2";
	mkdir(exemplo);

	//corpo0
	m[0]=0.3*mSol;
	printf("%lf\n",m[0]);

	//corpo1

	m[1]=mJupiter;
	double periodo=7.5;
	z[1*dimensao+0]=G*m[0]*pow(periodo,2)/(4*pow(PI,2));
	z[1*dimensao+0]=pow(z[1*dimensao+0],((double)1.0/3));
	z[4*dimensao+1]=m[1]*sqrt(G*m[0]/z[1*dimensao+0]);

	//corpo2
	m[2]=mTerra;

	//posicao inicial corpo 2:
	//a2=semi-eixo maior, b=semi-eixo menor, c= distância focal,
	double c=e*a2;
	double b=sqrt(pow(a2,2)-pow(c,2));
	z[2*dimensao+0]=a2*(1-e);

	double P=sqrt(4*pow(PI,2)*pow(a2,3)/(G*m[0]));
	double N=2*PI/P;
	z[5*dimensao+1]=m[2]*a2*b*N/z[2*dimensao+0];
}

void oito(void)
{
	//corpo0
	m[0]=1;
	z[0*dimensao+0]=0.9700436;
	z[0*dimensao+1]=-0.24308753;
	z[3*dimensao+0]=0.466203685;
	z[3*dimensao+1]=0.43236573;

	//corpo1
	m[1]=1;
	z[1*dimensao+0]=-0.9700436;
	z[1*dimensao+1]=+0.24308753;
	z[4*dimensao+0]=0.466203685;
	z[4*dimensao+1]=0.43236573;

	//corpo2
	m[2]=1;
	z[5*dimensao+0]=-2*0.466203685;
	z[5*dimensao+1]=-2*0.43236573;
}

void artigo1(void)
{
	//exemplo="/media/alan/Iomega HDD/artigo1-1";
	//exemplo="artigo1";
	exemplo="C:\\dados";
	mkdir(exemplo);
	exemplo="C:\\dados\\artigo1";
	mkdir(exemplo);

	//corpo0
	m[0]=mSol;

	//corpo1

	m[1]=mJupiter;
	double periodo=4.0;
	z[1*dimensao+0]=G*m[0]*pow(periodo,2)/(4*pow(PI,2));
	z[1*dimensao+0]=pow(z[1*dimensao+0],(double)(1.0/3));
	z[4*dimensao+1]=m[1]*sqrt(G*m[0]/z[1*dimensao+0]);


	//corpo2
	m[2]=mTerra;

	double c=e*a2;
	double b=sqrt(pow(a2,2)-pow(c,2));
	z[2*dimensao+0]=a2*(1-e);

	double P=sqrt(4*pow(PI,2)*pow(a2,3)/(G*m[0]));
	double N=2*PI/P;
	z[5*dimensao+1]=m[2]*a2*b*N/z[2*dimensao+0];
}



void coordenadasCentroMassa(void)
{
	int i,j;

	double M=0;
	double *rcm;
	rcm =(double *)calloc(dimensao, sizeof(double));

	for(i=0;i<corpos;i++)
	{
		M=M+m[i];
	}

	for(i=0;i<dimensao;i++)
	{
		rcm[i]=0;
		for(j=0;j<corpos;j++)
		{
			rcm[i]=rcm[i]+z[j*dimensao+i]*m[j];
		}
		rcm[i]=rcm[i]/M;
	}

	//transformacao
	for(i=0;i<corpos;i++)
	{
		for(j=0;j<dimensao;j++)
		{
			z[i*dimensao+j]= z[i*dimensao+j]-rcm[j];
		}
	}
}

void criarArquivos(void)
{
	char *raiz,*diretorio,*aux,*barra;
	int i;

	raiz=(char*)calloc(100,sizeof(char));
	aux=(char*)calloc(100,sizeof(char));
	diretorio=(char*)calloc(100,sizeof(char));

	if(!strcmp(SO,"w"))
	{
		barra="\\";
	}
	if(!strcmp(SO,"l"))
	{
		barra="/";
	}

	strcpy(raiz,exemplo);
	strcat(raiz,barra);
	//mkdir(raiz,S_IRWXU);
	mkdir(raiz);
	sprintf(aux,"%0.0lf",tFinal);
	strcat(raiz,aux);
	strcat(raiz,barra);
	//mkdir(raiz,S_IRWXU);
	mkdir(raiz);

	if(orbitas)
	{
		for(i=0;i<corpos;i++)
		{
			strcpy(diretorio,raiz);
			strcat(diretorio,"r");
			sprintf(aux,"%d",(i+1));
			strcat(diretorio,aux);
			strcat(diretorio,".txt");
			f[i] = fopen(diretorio,"w");
			if(f[i] == NULL)
			{
				printf("Erro ao abrir arquivo!!!\n");
				exit(1);
			}
			strcpy(diretorio,raiz);
			strcat(diretorio,"p");
			sprintf(aux,"%d",(i+1));
			strcat(diretorio,aux);
			strcat(diretorio,".txt");
			f[i+corpos] = fopen(diretorio,"w");
			if(f[i+corpos] == NULL)
			{
				printf("Erro ao abrir arquivo!!!\n");
				exit(1);
			}
		}
	}

	strcpy(diretorio,raiz);
	strcat(diretorio,"tempo");
	strcat(diretorio,".txt");
	f[2*corpos] = fopen(diretorio,"w");
	if(f[2*corpos] == NULL)
	{
		printf("Erro ao abrir arquivo5!!!\n");
		exit(1);
	}

	strcpy(diretorio,raiz);
	strcat(diretorio,"megno");
	strcat(diretorio,".txt");
	f[2*corpos+1] = fopen(diretorio,"w");
	if(f[2*corpos+1] == NULL)
	{
		printf("Erro ao abrir arquivo megno!!!\n");
		exit(1);
	}


	free(aux);
	free(raiz);
	free(diretorio);
}

void F3(double *dx, double *z, double *distancia, double *m, int dimensao)
{
	dx[0]=z[6]/m[0];
	dx[1]=z[7]/m[0];
	dx[2]=z[8]/m[1];
	dx[3]=z[9]/m[1];
	dx[4]=z[10]/m[2];
	dx[5]=z[11]/m[2];

	calcularDistancia(z,distancia,dimensao,corpos);

	dx[6]=      (G*m[0]*m[1]*(z[2]-z[0]))/(double)(pow(distancia[0*corpos+1],3));
	dx[6]=dx[6]+(G*m[0]*m[2]*(z[4]-z[0]))/(double)(pow(distancia[0*corpos+2],3));
	dx[7]=      (G*m[0]*m[1]*(z[3]-z[1]))/(double)(pow(distancia[0*corpos+1],3));
	dx[7]=dx[7]+(G*m[0]*m[2]*(z[5]-z[1]))/(double)(pow(distancia[0*corpos+2],3));

	dx[8]=      (G*m[1]*m[0]*(z[0]-z[2]))/(double)(pow(distancia[1*corpos+0],3));
	dx[8]=dx[8]+(G*m[1]*m[2]*(z[4]-z[2]))/(double)(pow(distancia[1*corpos+2],3));
	dx[9]=      (G*m[1]*m[0]*(z[1]-z[3]))/(double)(pow(distancia[1*corpos+0],3));
	dx[9]=dx[9]+(G*m[1]*m[2]*(z[5]-z[3]))/(double)(pow(distancia[1*corpos+2],3));

	dx[10]=       (G*m[2]*m[0]*(z[0]-z[4]))/(double)(pow(distancia[2*corpos+0],3));
	dx[10]=dx[10]+(G*m[2]*m[1]*(z[2]-z[4]))/(double)(pow(distancia[2*corpos+1],3));
	dx[11]=       (G*m[2]*m[0]*(z[1]-z[5]))/(double)(pow(distancia[2*corpos+0],3));
	dx[11]=dx[11]+(G*m[2]*m[1]*(z[3]-z[5]))/(double)(pow(distancia[2*corpos+1],3));

	jacobiana3C(J,m[0],m[1],m[2],z[0], z[1], z[2], z[3], z[4], z[5],tamanho2);

	produtoMatrizVetor(J,z+12,tamanho2,dx+tamanho2);

	produtoEscalar(dx+12,z+12,12,temp3);
	produtoEscalar(z+12,z+12,12,temp4);
	dx[24]=*temp3/(*temp4)*t;
	if(t==0)
	{
		t=1;
	}
	dx[25]=2*z[24]/t;

}

void calcularDistancia(double *z,double *distancia, int dimensao, int corpos)
{
	int i,j,k;
	double temp;

	for(i=0;i<corpos;i++)
	{
		for(j=0;j<corpos;j++)
		{
			distancia[i*corpos+j]=0;
			for(k=0;k<dimensao;k++)
			{
				temp=(z[i*dimensao+k]-z[j*dimensao+k]);
				temp=temp*temp;
				distancia[i*corpos+j]=distancia[i*corpos+j]+temp;
			}
			distancia[i*corpos+j]=(double)sqrt((double)distancia[i*corpos+j]);
		}
	}
}

void jacobiana3C(double* j, double m1,double m2,double m3,double r1x,double r1y,double r2x,double r2y,double r3x,double r3y,int tamanho)
{
	double G2=(double)G;
	double t70 = r1x-r2x;
	double t71 = r1y-r2y;
	double t72 = r1x-r3x;
	double t73 = r1y-r3y;
	double t78 = r1x*2;
	double t83 = r2x*2;
	double t74 = t78-t83;
	double t75 = pow(t70,2);
	double t76 = pow(t71,2);
	double t77 = t75+t76;
	double t85 = r3x*2;
	double t79 = -1*t78+t85;
	double t80 = pow(t72,2);
	double t81 = pow(t73,2);
	double t82 = t80+t81;
	double t84 = pow((1/t77),(5.0/2.0));
	double t86 = r1y*2;
	double t87 = pow((1/t82),(5.0/2.0));
	double t88 = pow((1/t77),(3.0/2.0));
	double t89 = pow(t74,2);
	double t90 = G2 * m1 * m2 * t84 * t89* 3.0/4.0;
	double t91 = r2y*2;
	double t92 = t86-t91;
	double t93 = G2*m1*m2*t74*t84*t92*3.0/4.0;
	double t94 = pow((1.0/t82),(3.0/2.0));
	double t95 = t78-t85;
	double t96 = r3y*2.0;
	double t97 = t86-t96;
	double t98 = G2*m1*m3*t87*t95*t97*(3.0/4.0);
	double t99 = t93+t98;
	double t100 = G2*m1*m2*t88;
	double t101 = G2*m1*m3*t94;
	double t102 = pow(t92,2);
	double t103 = G2*m1*m2*t102*t84*(3.0/4.0);
	double t104 = pow(t97,2);
	double t105 = G2*m1*m3*t104*t87*(3.0/4.0);
	double t106 = t100-t90;
	double t107 = r2x-r3x;
	double t108 = r2y-r3y;
	double t109 = t83-t85;
	double t110 = pow(t107,2);
	double t111 = pow(t108,2);
	double t112 = t110+t111;
	double t113 = pow((1.0/t112),(5.0/2.0));
	double t114 = pow((1.0/t112),(3.0/2.0));
	double t115 = pow(t109,2);
	double t116 = G2*m2*m3*t113*t115*(3.0/4.0);
	double t117 = t91-t96;
	double t118 = G2*m2*m3*t109*t113*t117*(3.0/4.0);
	double t119 = t100-t103;
	double t120 = t118+t93;
	double t121 = G2*m2*m3*t114;
	double t122 = pow(t117,2);
	double t123 = G2*m2*m3*t113*t122*(3.0/4.0);
	double t124 = pow(t95,2);
	double t127 = G2*m1*m3*t124*t87*(3.0/4.0);
	double t125 = t101-t127;
	double t126 = -1*t116+t121;
	double t128 = t101-t105;
	double t129 = t121-t123;
	double t130 = t118+t98;

	int colunas=3*2*2;

	j[0*colunas+6]=1/m1;
	j[1*colunas+7]=1/m1;
	j[2*colunas+8]=1/m2;
	j[3*colunas+9]=1/m2;
	j[4*colunas+10]=1/m3;
	j[5*colunas+11]=1/m3;
	j[6*colunas+0]=t90-G2*m1*m2*t88-G2*m1*m3*t94+G2*m1*m3*pow(t79,2)*t87*(3.0/4.0);
	j[6*colunas+1]=t99;
	j[6*colunas+2]=t106;
	j[6*colunas+3]=-t93;
	j[6*colunas+4]=t125;
	j[6*colunas+5]=-t98;
	j[7*colunas+0]=j[0*colunas+1];
	j[7*colunas+1]=-t100-t101+t103+t105;
	j[7*colunas+2]=-t93;
	j[7*colunas+3]=t119;
	j[7*colunas+4]=-t98;
	j[7*colunas+5]=t128;
	j[8*colunas+0]= t106;
	j[8*colunas+1]=-t93;
	j[8*colunas+2]=-t100+t116+t90-G2*m2*m3*t114;
	j[8*colunas+3]=t120;
	j[8*colunas+4]=t126;
	j[8*colunas+5]=-t118;
	j[9*colunas+0]=-t93;
	j[9*colunas+1]=t119;
	j[9*colunas+2]=t120;
	j[9*colunas+3]=-t100+t103-t121+t123;
	j[9*colunas+4]=-t118;
	j[9*colunas+5]=t129;
	j[10*colunas+0]=t125;
	j[10*colunas+1]=-t98;
	j[10*colunas+2]=t126;
	j[10*colunas+3]=-t118;
	j[10*colunas+4]=-t101+t116-t121+t127;
	j[10*colunas+5]=t130;
	j[11*colunas+0]=-t98;
	j[11*colunas+1]=t128;
	j[11*colunas+2]=-t118;
	j[11*colunas+3]=t129;
	j[11*colunas+4]=t130;
	j[11*colunas+5]=-t101+t105-t121+t123;
}

void produtoEscalar(double* v,double* v2,int tamanho, double *resultado)
{
	int i;

	resultado[0]=0;
	for(i=0;i<tamanho;i++)
	{
		resultado[0]+=v[i]*v2[i];
	}
}

void copiarVetor(double* v,double* v2)
{
	int i;

	for(i=0;i<tamanho;i++)
	{
		v[i]=v2[i];
	}
}

void somaVetores(double* v,double* v2,int tamanho, double *resultado)
{
	int i;

	for(i=0;i<tamanho;i++)
	{
		resultado[i]=v[i]+v2[i];
	}
}

void diferencaVetores(double* v,double* v2,int tamanho, double *resultado)
{
	int i;

	for(i=0;i<tamanho;i++)
	{
		resultado[i]=v[i]-v2[i];
	}
}

void produtoEscalarVetor(double* v,double a,int tamanho,double *resultado)
{
	int i;

	for(i=0;i<tamanho;i++)
	{
		resultado[i]=a*v[i];
	}
}

void produtoMatrizVetor(double *J,double* delta,int tamanho, double *resultado)
{
	int i,j;
	double temp;

	for(i=0;i<tamanho;i++)
	{
		temp=J[i*tamanho]*delta[0];
		for(j=1;j<tamanho;j++)
		{
			temp+=J[i*tamanho+j]*delta[j];
		}
		resultado[i]=temp;
	}
}

void escreverArquivos(void)
{
	int j,k;
	char aux[100];

	if(orbitas)
	{
		for(k=0;k<corpos;k++)
		{
			for(j=0;j<dimensao;j++)
			{
				//arquivo txt
				//posicao
				sprintf(aux,"%0.20lf ",z[k*dimensao+j]);
				fputs(aux,f[k]);
				//momento
				sprintf(aux,"%0.20lf ",z[(k+corpos)*dimensao+j]);
				fputs(aux,f[k+corpos]);
			}
			fputs("\n",f[k]);
			fputs("\n",f[k+corpos]);
		}
	}
	sprintf(aux,"%0.15lf ",z[tamanho-1]/t);
	fputs(aux,f[2*corpos+1]);
	fputs("\n",f[2*corpos+1]);
}

void imprimir(double* v,int tamanho)
{
	int i;
	for(i=0;i<tamanho;i++)
	{
		printf("%0.10lf\n",v[i]);
	}
	puts("");
}

double calculaNorma(double *v,int dimensao)
{
	double temp=0;
	int i;

	for(i=0;i<dimensao;i++)
	{
		temp+=v[i]*v[i];
	}
	temp=sqrt(temp);
	return temp;
}

void criarMatrizRK45(void)
{
	int stages=6;
	a=(double*)calloc(stages*stages, sizeof(double));
	b4=(double*)calloc(stages, sizeof(double));
	b5=(double*)calloc(stages, sizeof(double));
	c2=(double*)calloc(stages, sizeof(double));

	a[1*stages+0]=1.0/4;
	a[2*stages+0]=3.0/32;
	a[2*stages+1]=9.0/32;
	a[3*stages+0]=1932.0/2197;
	a[3*stages+1]=-7200.0/2197;
	a[3*stages+2]=7296.0/2197;
	a[4*stages+0]=439.0/216;
	a[4*stages+1]=-8.0;
	a[4*stages+2]=3680.0/513;
	a[4*stages+3]=-845.0/4104;
	a[5*stages+0]=-8.0/27;
	a[5*stages+1]=2.0;
	a[5*stages+2]=-3544.0/2565;
	a[5*stages+3]=1859.0/4104;
	a[5*stages+4]=-11.0/40;

	b4[0]=25.0/216;
	b4[2]=1408.0/2565;
	b4[3]=2197.0/4104;
	b4[4]=-1.0/5;

	b5[0]=16.0/135;
	b5[2]=6656.0/12825;
	b5[3]=28561.0/56430;
	b5[4]=-9.0/50;
	b5[5]=2.0/55;

	//independente
	c2[0]=0.0;
	c2[1]=1.0/4;
	c2[2]=3.0/8;
	c2[3]=12.0/13;
	c2[4]=1.0;
	c2[5]=1.0/2;
}

void finalizar(void)
{
	if(orbitas)
	{
		fclose(f[0]);
		fclose(f[1]);
		fclose(f[2]);
		fclose(f[3]);
		fclose(f[4]);
		fclose(f[5]);
	}
	fclose(f[6]);
	fclose(f[7]);

	free(distancia);
	free(m);
	free(a);
	free(b4);
	free(b5);
	free(c2);
	free(K1);
	free(K2);
	free(K3);
	free(K4);
	free(K5);
	free(K6);
	free(J);
	free(temp);
	free(norma);
	free(Y);
	free(y);
	free(z);
}
